#ifndef _COUPLING_COEFF_
#define _COUPLING_COEFF_
#include "Wigner_signs.h"
//#include "gsl/gsl_sf_coupling.h"
#include <cmath>
#include <iostream>
#include <stdlib.h>
#include <vector>

using namespace Wigner_signs;
//
// extern struct {
//   int n;
//   double f;
//   long i;
// } fact_table_coupling[];
// extern struct {
//   int n;
//   double f;
//   long i;
// } doub_fact_table_coupling[];
//~~~~~~~~start cg coeff ~~~~~~~~~~~~~~~~~~~//
struct j2m2_state {
  // jama jbmb -> JM
  int j2a;
  int m2a;
  int size;
  // int index;
};

struct jajbJ_state {
  // ja jb -> J
  int j2a;
  int j2b;
  int J2;
  int index; // index of ja jb -> J state
};
struct TrK_state {
  // K is the index of jajbJ state
  // the state of k1 k2 k3
  int k1;
  int k2;
  int k3;
  // int index;
};

class J2M2 {
public:
  std::vector<std::vector<int>> vec; // vec[j2a][m2a+j2a] num of m = j2a+1
  std::vector<j2m2_state> state;
  //#pragma acc routine seq
  void init(int &j2a_max);
};

class TrK_cg {
  // the state of k1 k2 k3
  // index of TrK =TrK.vec[k1][exindex2(k2)][exindex3(k3)]
public:
  std::vector<std::vector<std::vector<int>>> vec;
  std::vector<TrK_state> state;
  std::vector<std::vector<std::vector<int>>> exindex3; // exchange index for k3
  std::vector<std::vector<std::vector<int>>> inv_exindex3;
  //#pragma acc routine seq
  void init(J2M2 &A);
};

class CG_coeff {
public:
  int init_info = 1;
  CG_coeff();
  //#pragma acc routine seq
  double cal(int j2a, int j2b, int J2, int m2a, int m2b, int M2) const;
  int j2a_max;
  // int init_info = 0;
  std::vector<double> vec; // cg = vec(index)
  J2M2 J2M2_t;
  TrK_cg TrKcg_t;
  void init_cal(J2M2 &A, TrK_cg &B);
  void init(int j2a_max) {
    this->j2a_max = j2a_max;
    J2M2_t.init(j2a_max);
    TrKcg_t.init(J2M2_t);
    // std::cout << "!! @@ !!" << TrKcg_t.state.size() << std::endl;
    this->init_cal(this->J2M2_t, this->TrKcg_t);
    this->init_info = 1;
  }
  // quick find cg value
  //#pragma acc routine seq
  double quick(int j2a, int j2b, int J2, int m2a, int m2b, int M2) const;
};

//~~~~~~~~end cg coeff ~~~~~~~~~~~~~~~~~~~//

// ~~~~~~~~for six J coeff~~~~~~~~~~~~~~~~~~//

class JaJbJ {
  // coupling state of <jajb|J>
public:
  std::vector<std::vector<std::vector<int>>> vec;
  std::vector<jajbJ_state> state;
  int Tri(int &j2a, int &j2b,
          int &j2c); // triadic condition for ja,jb,jc: -1 or 1
                     //#pragma acc routine seq
  void init(int &J2a_max);
};
class TrK {
  // the state of k1 k2 k3
  // index of TrK =TrK.vec[k1][exindex2(k2)][exindex3(k3)]
public:
  std::vector<std::vector<std::vector<int>>> vec;
  std::vector<TrK_state> state;
  std::vector<std::vector<int>> exindex2;              // exchange index for k2
  std::vector<std::vector<std::vector<int>>> exindex3; // exchange index for k3
  // std::vector<std::vector<int>> inv_exindex2;
  // std::vector<std::vector<std::vector<int>>> inv_exindex3;
  void init(JaJbJ &A);
};

class SixJ_coeff {
public:
  SixJ_coeff();
  // double cal(int j2a, int j2b, int J2ab, int j2c, int J2, int J2bc) const {
  //   return gsl_sf_coupling_6j(j2a, j2b, J2ab, j2c, J2, J2bc);
  // }
  int init_info = 1;
#pragma acc routine seq
  double cal(int j2a, int j2b, int J2ab, int j2c, int J2, int J2bc) const;
  // for quick calculation
  // store 6j data

  int j2a_max;
  // int init_info;
  std::vector<double> vec; // 6j = sixJ_coeff.vec(TrK.index)
  JaJbJ JaJbJ_t;
  TrK TrK_t;
  void init_cal(JaJbJ &A, TrK &B);
  void init(int j2a_max) {
    this->j2a_max = j2a_max;
    JaJbJ_t.init(j2a_max);
    TrK_t.init(this->JaJbJ_t);
    this->init_cal(this->JaJbJ_t, this->TrK_t);
    this->init_info = 1;
  }
  // quick find 6j value
#pragma acc routine seq
  double quick(int j2a, int j2b, int J2ab, int j2c, int J2, int J2bc) const;
};

// ~~~~~~~~for nine J coeff~~~~~~~~~~~~~~~~~~//
class TrK_ninJ {
public:
  std::vector<std::vector<std::vector<int>>> vec;
  std::vector<TrK_state> state;
  std::vector<std::vector<std::vector<int>>> exindex3;
  void init(JaJbJ &A);
};
class NinJ_coeff {
public:
  NinJ_coeff();
  NinJ_coeff(SixJ_coeff &sixJ_t) : sixJ(sixJ_t){};
  // double cal(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd, int
  // J2ac,
  //            int J2bd, int J2) const {
  //   return gsl_sf_coupling_9j(j2a, j2b, J2ab, j2c, j2d, J2cd, J2ac, J2bd,
  //   J2);
  // };
  int init_info = 1;
#pragma acc routine seq
  double cal(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd, int J2ac,
             int J2bd, int J2) const;
  int j2a_max;
  // int init_info;
  SixJ_coeff sixJ;
  std::vector<double> vec; // 9j = ninJ_coeff.vec(TrK.index)
  JaJbJ JaJbJ_t;
  TrK_ninJ TrKninJ_t;
  void init_cal(JaJbJ &A, TrK_ninJ &B);
  void init(int j2a_max, SixJ_coeff &sixJ_t) {
    this->j2a_max = j2a_max;
    JaJbJ_t.init(j2a_max);
    TrKninJ_t.init(this->JaJbJ_t);
    this->init_cal(this->JaJbJ_t, this->TrKninJ_t);
    this->init_info = 1;
    this->sixJ = sixJ_t;
  }
#pragma acc routine seq
  double quick(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd, int J2ac,
               int J2bd, int J2) const;
#pragma acc routine seq
  double quick(int &j2a, int &j2b, int &J2ab, int &j2c, int &j2d, int &J2cd,
               int &J2ac, int &J2bd, int &J2) const;
  // double quick2(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd,
  //               int J2ac, int J2bd, int J2) const;

private:
  int phase(int &n) const { return n % 2 ? -1 : 1; }
};

// ~~~~~~~~end nine J coeff~~~~~~~~~~~~~~~~~~//
//~~~~~~~~~ begin  binomial coeff ~~~~~~~~~~~//
class Binom {
public:
  int n_max;
  int init_info;
  std::vector<std::vector<std::vector<double>>>
      vec;                     // vec[n][k1][k2] n!!/(k1!!*k2!!)
  std::vector<double> fac2vec; // n!! = fac2vec[n]
  std::vector<double> facvec;  // n! = facvec[n]

  Binom() {
    n_max = 0;
    init_info = -1;
    // vec.clear();
    // fac2vec.clear();
    // facvec.clear();
    //#pragma acc enter data copyin(this)
    //#pragma acc enter data copyin(n_max)
    //#pragma acc enter data copyin(vec)
    //#pragma acc enter data copyin(facvec)
    //#pragma acc enter data copyin(fac2vec)
  }
  //   ~Binom() {
  // //#pragma acc exit data delete (n_max, init_info)
  // //#pragma acc exit data delete (this)
  //     vec.clear();
  //     fac2vec.clear();
  //     facvec.clear();
  //   }
  //#pragma acc routine seq
  void fac2veccal(int &n);
  //#pragma acc routine seq
  void facveccal(int &n);
#pragma acc routine seq
  double fac_cal(int n);
#pragma acc routine seq
  double fac2_cal(int n);
#pragma acc routine seq
  double binom_cal(int n, int k1, int k2);
  //#pragma acc routine seq
  double init_cal(int &n_max);
  //#pragma acc routine seq
  void init(int &n_max){
      // this->fac2veccal(n_max);
      // this->facveccal(n_max);
      // this->init_cal(n_max);
      // this->init_info = 1;
      //#pragma acc enter data copyin(this)
      //#pragma acc enter data copyin(n_max)
      //#pragma acc enter data copyin(vec)
      //#pragma acc enter data copyin(facvec)
      //#pragma acc enter data copyin(fac2vec)
  };
  void update_device() {
    //#pragma acc enter data copyin(n_max)
    //#pragma acc enter data copyin(vec)
    //#pragma acc enter data copyin(facvec)
    //#pragma acc enter data copyin(fac2vec)
  }
  void update_host() {
    //#pragma acc enter data copyin(n_max)
    //#pragma acc enter data copyin(vec)
    //#pragma acc enter data copyin(facvec)
    //#pragma acc enter data copyin(fac2vec)
  }
  //#pragma acc routine seq
  double quick(int &i, int &j, int &k) { return this->vec[i][j][k]; }
  //#pragma acc routine seq
  double quick(int i, int j, int k) { return this->vec[i][j][k]; }
};
//~~~~~~~~~ end  binomial coeff ~~~~~~~~~~~//

//~~~~~~~~~ begin  harmonic-oscillator brackets ~~~~~~~~~~~//
class Hobra {
public:
  Hobra(Binom &binom_t, CG_coeff &cg_t, NinJ_coeff &ninJ_t)
      : binom_p(binom_t), cg_p(cg_t), ninJ_p(ninJ_t){};
  Binom &binom_p;
  NinJ_coeff &ninJ_p;
  CG_coeff &cg_p;

  std::vector<std::vector<std::vector<double>>>
      vec;                     // vec[n][k1][k2] n!!/(k1!!*k2!!)
  std::vector<double> fac2vec; // n!! = fac2vec[n]
  std::vector<double> facvec;  // n! = facvec[n]

  void fac2veccal(int n);
  void facveccal(int n);
  double binom_cal(int n, int k1, int k2) const;
  void init_binom(int n_max_fac, int n_max_fac_2);

  void init(Binom &binom_t, CG_coeff &cg_t, NinJ_coeff &ninJ_t) {
    this->binom_p = binom_t;
    this->cg_p = cg_t;
    this->ninJ_p = ninJ_t;
    // this->init_binom(170, 297);
  }
  void init_2(Binom &binom_t, CG_coeff &cg_t, NinJ_coeff &ninJ_t) {
    this->binom_p = binom_t;
    this->cg_p = cg_t;
    this->ninJ_p = ninJ_t;
    this->init_binom(170, 297);
  }

  int init_info = this->cg_p.init_info * this->ninJ_p.init_info;
#pragma acc routine seq
  double fac_cal(int n);
#pragma acc routine seq
  double fac2_cal(int n);
#pragma acc routine seq
  double cal(int &E, int &L, int &e, int &l, int &e1, int &l1, int &e2, int &l2,
             int &lambda, double &d) const;

#pragma acc routine seq
  double cal(int E, int L, int e, int l, int e1, int l1, int e2, int l2,
             int lambda, double d, int k) const {
    return this->cal(E, L, e, l, e1, l1, e2, l2, lambda, d);
  }

#pragma acc routine seq
  double cal2(int &E, int &L, int &e, int &l, int &e1, int &l1, int &e2,
              int &l2, int &lambda, double &d) const;

  void test() {
    std::cout << " init info = " << this->binom_p.init_info << std::endl;
  }

private:
#pragma acc routine seq
  double G(int &e1, int &l1, int &e_a, int &l_a, int &e_b, int &l_b) const;
  double G_2(int &e1, int &l1, int &e_a, int &l_a, int &e_b, int &l_b) const;
  //#pragma acc routine seq
  int phase_cal(int n) const { return n % 2 ? -1 : 1; }
};
//~~~~~~~~~ end  harmonic-oscillator brackets ~~~~~~~~~~~//

//~~~~~~~~~ some functions~~~~~~~~~~~~~~~~~~//
class Pre_coeff {
public:
  // Pre_coeff(){};
  // bino
  std::vector<std::vector<std::vector<double>>>
      vec_bino;                // vec[n][k1][k2] n!!/(k1!!*k2!!)
  std::vector<double> fac2vec; // n!! = fac2vec[n]
  std::vector<double> facvec;  // n! = facvec[n]
  void init_Binom();
#pragma acc routine seq
  double fac_cal(int n);
#pragma acc routine seq
  double fac_cal_this(int n);
#pragma acc routine seq
  double fac2_cal(int n);
#pragma acc routine seq
  double fac2_cal_this(int n);
#pragma acc routine seq
  double binom_cal_this(int n, int k1, int k2);
#pragma acc routine seq
  double binom_cal(int n, int k1, int k2) {
    return vec_bino[n][k1][k2];
  }
  //~CG
  J2M2 J2M2_t;
  int j2a_max_cg;
  std::vector<double> vec_cg;
  TrK_cg TrKcg_t;
  void get_CG(CG_coeff cg_t) {
    this->j2a_max_cg = cg_t.j2a_max;
    this->J2M2_t = cg_t.J2M2_t;
    this->vec_cg = cg_t.vec;
    this->TrKcg_t = cg_t.TrKcg_t;
  }
#pragma acc routine seq
  double cg_quick(int j2a, int j2b, int J2, int m2a, int m2b, int M2);
#pragma acc routine seq
  double cg_cal_this(int j2a, int j2b, int J2, int m2a, int m2b, int M2);
  // ~~~~~~~~~~~~~~~~~~~ 6j ~~~~~~~~~~~~~~~~~~~ //
  std::vector<double> vec_6j; // 6j = sixJ_coeff.vec(TrK.index)
  JaJbJ JaJbJ_t;
  TrK TrK_t;
  int j2a_max_6j;
  void get_SixJ(SixJ_coeff sixJ_t) {
    this->j2a_max_6j = sixJ_t.j2a_max;
    this->JaJbJ_t = sixJ_t.JaJbJ_t;
    this->TrK_t = sixJ_t.TrK_t;
    this->vec_6j = sixJ_t.vec;
  }
#pragma acc routine seq
  double SixJ_quick(int j2a, int j2b, int J2ab, int j2c, int J2, int J2bc);
#pragma acc routine seq
  double sixJ_this(int j2a, int j2b, int J2ab, int j2c, int J2, int J2bc);
  // ~~~~~~~~~~~~~~~~~~~ 9j ~~~~~~~~~~~~~~~~~~~ //
  // int j2a_max_9j;
  // std::vector<double> vec_9j; // 9j = ninJ_coeff.vec(TrK.index)
  // JaJbJ JaJbJ_t_9j;
  // TrK_ninJ TrKninJ_t;
  // void get_NinJ(NinJ_coeff ninJ_t) {
  //   this->j2a_max_9j = ninJ_t.j2a_max;
  //   this->JaJbJ_t_9j = ninJ_t.JaJbJ_t;
  //   this->TrKninJ_t = ninJ_t.TrKninJ_t;
  //   this->vec_9j = ninJ_t.vec;
  // }
#pragma acc routine seq
  double NinJ_quick(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd,
                    int J2ac, int J2bd, int J2);
//
// ~~~~~~~~~~~~~~~~~~~ Hobra ~~~~~~~~~~~~~~~~~~~ //
#pragma acc routine seq
  double hobra_cal(int &E, int &L, int &e, int &l, int &e1, int &l1, int &e2,
                   int &l2, int &lambda, double &d);

private:
#pragma acc routine seq
  double G(int &e1, int &l1, int &e_a, int &l_a, int &e_b, int &l_b);
  int phase(int n) const { return n % 2 ? -1 : 1; }
};
/// gauss-legendre quadrature
// void gauleg(const double x1, const double x2, vector<double> &x,
// vector<double> &w);

#endif
